home *** CD-ROM | disk | FTP | other *** search
/ HPAVC / HPAVC CD-ROM.iso / pc / SNNSV32.ZIP / SNNSv3.2 / kernel / sources / matrix.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-04-25  |  17.5 KB  |  683 lines
  1. /*****************************************************************************
  2.   FILE           : matrix.c
  3.   SHORTNAME      : 
  4.   SNNS VERSION   : 3.2
  5.  
  6.   PURPOSE        : SNNS-Kernel standard matrix operations
  7.   NOTES          :
  8.  
  9.   AUTHOR         : Michael Vogt
  10.   DATE           : 10.02.92
  11.  
  12.   CHANGED BY     : Sven Doering
  13.   IDENTIFICATION : @(#)matrix.c    1.10 3/15/94
  14.   SCCS VERSION   : 1.10
  15.   LAST CHANGE    : 3/15/94
  16.  
  17.              Copyright (c) 1990-1994  SNNS Group, IPVR, Univ. Stuttgart, FRG
  18.  
  19. ******************************************************************************/
  20.  
  21. /************************************************************************/
  22. /* included headers:                            */
  23. /************************************************************************/
  24.  
  25. #include <malloc.h>
  26. #include <stdio.h>
  27. #include <math.h>
  28. #include "matrix.ph"
  29. #include "glob_typ.h"
  30.  
  31.  
  32. /************************************************************************/
  33. /* functions:                                */
  34. /************************************************************************/
  35.  
  36. /************************************************************************/
  37. /* function RbfAllocMatrix:                        */
  38. /* allocate matrix m with r rows and c columns                */
  39. /* returns 0 if impossible, 1 otherwise                    */
  40. /************************************************************************/
  41.  
  42. int    RbfAllocMatrix(int r, int c, RbfFloatMatrix *m)
  43. {
  44.     int    rc;            /* return code            */
  45.     int    i;            /* index variable        */
  46.  
  47.     m -> field = (float *) malloc(r * c * sizeof(float));
  48.     m -> r_pt = (float **) malloc(r * sizeof(float *));
  49.     if (m -> field == (float *) NULL ||
  50.         m -> r_pt == (float **) NULL)
  51.     {
  52.         /* malloc was impossible, return 0 (reports error)    */
  53.         rc = 0;
  54.     }
  55.     else
  56.     {
  57.         /* malloc successfull, initialize matrix:        */
  58.         rc = 1;
  59.         m -> rows = r;
  60.         m -> columns = c;
  61.         for (i = 0; i<r; i++)
  62.             (m -> r_pt)[i] = &(m -> field)[i * c];
  63.     }
  64.  
  65.     return rc;            /* report status        */
  66. }
  67.  
  68. /************************************************************************/
  69. /* function RbfFreeMatrix:                        */
  70. /* deallocate matrix m                            */
  71. /************************************************************************/
  72.  
  73. void    RbfFreeMatrix(RbfFloatMatrix *m)
  74. {
  75.     free(m -> field);
  76.     free(m -> r_pt);
  77.     m -> rows = 0;
  78.     m -> columns = 0;
  79. }
  80.  
  81. /************************************************************************/
  82. /* function RbfClearMatrix:                        */
  83. /* set all elements of matrix m to value c                */
  84. /************************************************************************/
  85.  
  86. void RbfClearMatrix(RbfFloatMatrix *m, double c)
  87. {
  88.     int    count;
  89.     float    *ptoelement;
  90.  
  91.     ptoelement = m -> field;
  92.     count = (m -> rows)*(m -> columns);
  93.  
  94.     while(count--)
  95.     {
  96.         *ptoelement++ = c;
  97.     }
  98. }
  99.  
  100. /************************************************************************/
  101. /* function RbfSquareOfNorm:                        */
  102. /*                                                      */
  103. /************************************************************************/
  104.  
  105. float RbfSquareOfNorm(RbfFloatMatrix *m)
  106. {
  107.   int    i, j;
  108.   float norm = 0.0;
  109.  
  110.   for (i = m->rows -1 ; i>=0; i--)
  111.     {
  112.       for (j = m->columns -1 ; j>=0; j--)
  113.     norm += RbfMatrixGetValue(m, i ,j )*RbfMatrixGetValue(m, i, j);
  114.     };
  115.   return norm ;
  116. }
  117.  
  118. /************************************************************************/
  119. /* function RbfIdempotentMatrix:                                    */
  120. /* CAUTION: m must be a square matrix                            */
  121. /* set all elements of matrix m to value 0, the diagonal to 1          */
  122. /************************************************************************/
  123.  
  124. void RbfIdempotentMatrix(RbfFloatMatrix *m)
  125. {
  126. register int     i,j;
  127.  
  128. for (i = m->rows -1 ; i>=0; i--)
  129.   {
  130.     for (j = m->columns -1 ; j>=0; j--)
  131.       RbfMatrixSetValue(m,  i, j, 0);
  132.     RbfMatrixSetValue(m, i, i, 1);
  133.   };
  134.  
  135. }
  136.  
  137. /************************************************************************/
  138. /* function RbfMulScalarMatrix:                                    */
  139. /* set all elements of matrix m to their value multiplied with a          */
  140. /************************************************************************/
  141.  
  142. void RbfMulScalarMatrix(RbfFloatMatrix *m, float a)
  143. {
  144. register int     i,j;
  145.  
  146. for (i = m->rows -1 ; i>=0; i--)
  147.   {
  148.     for (j = m->columns -1 ; j>=0; j--)
  149.       RbfMatrixSetValue(m, i, j, a * RbfMatrixGetValue(m, i, j) );
  150.   };
  151.  
  152. }
  153.  
  154. /************************************************************************/
  155. /* function RbfSetMatrix:                        */
  156. /* set m1 to m2 (m1 = m2). m1 must be allocated with the same        */
  157. /* dimensions as m2                            */
  158. /************************************************************************/
  159.  
  160. void RbfSetMatrix(RbfFloatMatrix *m1, RbfFloatMatrix *m2)
  161. {
  162.     int    count;
  163.     float    *ptofrom;
  164.     float    *ptoto;
  165.  
  166. #ifdef    DEBUG_MODE
  167.     if (m1 -> rows != m2 -> rows ||
  168.         m1 -> columns != m2 -> columns)
  169.     {
  170.         ErrMess("RbfSetMatrix: incompatible matrixes.\n");
  171.     }
  172. #endif
  173.  
  174.     count = (m2 -> rows)*(m2 -> columns);
  175.     ptofrom = m2 -> field;
  176.     ptoto = m1 -> field;
  177.  
  178.     while(count--)
  179.     {
  180.         *ptoto++ = *ptofrom++;
  181.     }
  182. }
  183.  
  184. /************************************************************************/
  185. /* function RbfTranspMatrix:                        */
  186. /* set m1 to m2 transposed (m1 = m2T)                    */
  187. /* number of rows of m1 must be equal to number of columns of m2.    */
  188. /* number of columns of m1 must be equal to number of rows of m2.    */
  189. /************************************************************************/
  190.  
  191. void    RbfTranspMatrix(RbfFloatMatrix *m1, RbfFloatMatrix *m2)
  192. {
  193.     
  194.     int        r;
  195.     int        c;
  196.  
  197. #ifdef    DEBUG_MODE
  198.     if (m1 -> rows != m2 -> columns ||
  199.         m1 -> columns != m2 -> rows)
  200.         ErrMess("RbfTranspMatrix: incompatible matrixes.\n"); 
  201. #endif
  202.  
  203.     /* make m1 to become m2 transposed                */
  204.  
  205.     for (r = 0; r < m2 -> rows; r++)
  206.         for (c = 0; c < m2 -> columns; c++)
  207.         RbfMatrixSetValue(m1, c, r, RbfMatrixGetValue(m2, r, c));
  208. }
  209.  
  210. #ifdef RBF_MATRIX_LUDCOMP
  211.  
  212. /************************************************************************/
  213. /* function RbfLUDcmp:                            */
  214. /* builds the LU Decomposition of matrix m and stores it into m     */
  215. /* the permutation of rows is stored in vector <index> to use with    */
  216. /* the function RbfLUBksb for solving linear equations.            */
  217. /* (Algorithm taken from 'Numerical Recipts')                */
  218. /* returns 1 if succesfull, 0 if singular matrix, negative if kernel    */
  219. /* error                                */
  220. /************************************************************************/
  221.  
  222. static int RbfLUDcmp(RbfFloatMatrix *m, int *indx)
  223. {
  224.     register float        sum;
  225.     register float        dum;
  226.     register float        big;
  227.     register int        i, j, k, imax;
  228.     register float        temp;
  229.     register float        *vv;
  230.  
  231.     if ((vv = (float *) malloc(m -> rows * sizeof(float))) == NULL)
  232.     {
  233.         ErrMess("RbfLUDcmp: impossible to allocate helpmatrix.\n");
  234.         return KRERR_INSUFFICIENT_MEM;
  235.     }
  236.  
  237.     for (i = 0; i < m -> rows; i++)
  238.     {
  239.         big = 0.0;
  240.         for (j = 0; j < m -> rows; j++)
  241.         {
  242.         if ((temp = fabs(RbfMatrixGetValue(m, i, j))) > big)
  243.             big = temp;
  244.         }
  245.         if (big == 0.0)
  246.         {
  247.         free(vv);
  248.         return 0;
  249.         }
  250.         vv[i] = 1.0/big;
  251.     }
  252.     for (j = 0; j < m -> rows; j++)
  253.     {
  254.         for (i = 0; i < j; i++)
  255.         {
  256.         sum = RbfMatrixGetValue(m, i, j);
  257.         for (k = 0; k < i; k++)
  258.             sum -= RbfMatrixGetValue(m, i, k) *
  259.                RbfMatrixGetValue(m, k, j);
  260.         RbfMatrixSetValue(m, i, j, sum);
  261.         }
  262.         big = 0.0;
  263.         for (i = j; i < m -> rows; i++)
  264.         {
  265.         sum = RbfMatrixGetValue(m, i, j);
  266.         for (k = 0; k < j; k++)
  267.             sum -= RbfMatrixGetValue(m, i, k) *
  268.                RbfMatrixGetValue(m, k, j);
  269.         RbfMatrixSetValue(m, i, j, sum);
  270.         if ((dum = vv[i] * fabs(sum)) >= big)
  271.         {
  272.             big = dum;
  273.             imax = i;
  274.         }
  275.         }
  276.         if (j != imax)
  277.         {
  278.         for (k = 0; k < m -> rows; k++)
  279.         {
  280.             dum = RbfMatrixGetValue(m, imax, k);
  281.             RbfMatrixSetValue(m, imax, k, 
  282.             RbfMatrixGetValue(m, j, k));
  283.             RbfMatrixSetValue(m, j, k, dum);
  284.         }
  285.         dum = vv[imax];
  286.         vv[imax] = vv[j];
  287.         vv[j] = dum;
  288.         }
  289.         indx[j] = imax;
  290.         if (RbfMatrixGetValue(m, j, j) == 0.0)
  291.         {
  292.         fprintf(stderr,"RbfLUDcmp: seems to be a singular matrix\n");
  293.         free(vv);
  294.         return 0;
  295.         }
  296.         if (j != (m -> rows - 1))
  297.         {
  298.         dum = 1.0/RbfMatrixGetValue(m, j, j);
  299.         for (i = j+1; i < m -> rows; i++)
  300.             RbfMatrixSetValue(m, i, j,
  301.             RbfMatrixGetValue(m, i, j) * dum);
  302.         }
  303.     }
  304.     free(vv);
  305.     return 1;
  306. }
  307.  
  308. /************************************************************************/
  309. /* function RbfLUBksb                            */
  310. /* solves the linear equation LU*x = b, where LU is placed in m using   */
  311. /* the algorithm RbfLUDcmp. b is replaced by the solution x.        */
  312. /************************************************************************/
  313.  
  314. static void RbfLUBksb(RbfFloatMatrix *m, int *indx, float *b)
  315. {
  316.     register float    sum;
  317.     register int    i, ii=0, ip, j;
  318.  
  319.     for (i = 0; i < m -> rows; i++)
  320.     {
  321.         ip = indx[i];
  322.         sum = b[ip];
  323.         b[ip] = b[i];
  324.         if (ii)
  325.         {
  326.         for (j = ii-1; j < i; j++)
  327.             sum -= RbfMatrixGetValue(m, i, j) * b[j];
  328.         }
  329.         else if (sum != 0.0)
  330.         ii = i+1;
  331.         b[i] = sum;
  332.     }
  333.     for (i = (m -> rows - 1); i >= 0; i--)
  334.     {
  335.         sum = b[i];
  336.         for (j = i+1; j < m -> rows; j++)
  337.         sum -= RbfMatrixGetValue(m, i, j) * b[j];
  338.         b[i] = sum / RbfMatrixGetValue(m, i, i);
  339.     }
  340.     
  341. }
  342.  
  343. /************************************************************************/
  344. /* function RbfInvMatrix:                        */
  345. /* set m to inverse of m (m = m^-1).                    */
  346. /* returns 0 if impossible, 1 otherwise, negative if kernel error    */
  347. /* This function makes use of the Gauss-Jordan algorithm, which devides */
  348. /* the matrix m into two triangular matrixes m = l*r. m needs to be not */
  349. /* singulary! (Algorithm taken from "Numerical Recipts")        */
  350. /************************************************************************/
  351.  
  352. int    RbfInvMatrix(RbfFloatMatrix *m)
  353. {
  354.     register     int    i, j;
  355.     register    int    *indx;
  356.     register    float    *b;
  357.     register    int    tmp_err;
  358.     RbfFloatMatrix    help;
  359.  
  360. #ifdef    DEBUG_MODE
  361.     if (m -> rows != m -> columns)
  362.         ErrMess("RbfInvMatrix: matrix not of form N x N.\n");
  363. #endif
  364.  
  365.     if (!RbfAllocMatrix(m -> rows, m -> rows, &help) ||
  366.         (indx = (int *) malloc(m -> rows * sizeof(int))) == NULL ||
  367.         (b = (float *) malloc(m -> rows * sizeof(float))) == NULL)
  368.     {
  369.         ErrMess("RbfInvMatrix: impossible to allocate helpmatrix.\n");
  370.         return KRERR_INSUFFICIENT_MEM;
  371.     }
  372.  
  373.     RbfSetMatrix(&help, m);
  374.     if ((tmp_err = RbfLUDcmp(&help, indx)) != 1)
  375.     {
  376.         free(b);
  377.         free(indx);
  378.         RbfFreeMatrix(&help);
  379.         return tmp_err;
  380.     }
  381.  
  382.     for (j = 0; j < m -> rows; j++)
  383.     {
  384.         for (i = 0; i < m -> rows; i++)
  385.         b[i] = 0.0;
  386.         b[j] = 1.0;
  387.         RbfLUBksb(&help, indx, b);
  388.         for (i = 0; i < m -> rows; i++)
  389.         RbfMatrixSetValue(m, i, j, b[i]);
  390.     }
  391.  
  392.     free(b);
  393.     free(indx);
  394.     RbfFreeMatrix(&help);
  395.  
  396.     return 1;
  397. }
  398.  
  399. #else
  400.  
  401. /************************************************************************/
  402. /* function RbfInvMatrix:                        */
  403. /* set m to inverse of m (m = m^-1).                    */
  404. /* returns 0 if impossible, 1 otherwise, negative if kernel error    */
  405. /* This function makes use of the Gauss-Jordan algorithm, which devides */
  406. /* the matrix m into two triangular matrixes m = l*r. m needs to be not */
  407. /* singulary! (Algorithm taken from "Stoer: Numerische Mathematik 1")    */
  408. /************************************************************************/
  409.  
  410. int    RbfInvMatrix(m)
  411. RbfFloatMatrix    *m;
  412. {
  413.     RbfFloatMatrix    help;
  414.     int        i, j, r, k, n;
  415.     float        max, hr, hi;
  416.  
  417. #ifdef    DEBUG_MODE
  418.     if (m -> rows != m -> columns)
  419.         ErrMess("RbfInvMatrix: matrix not of form N x N.\n");
  420. #endif
  421.  
  422.     n = m -> rows;            /* n = dimension of matrix    */
  423.  
  424.     /* allocate helpmatrix: 0. row: field p[N] of algorithm        */
  425.     /*            1. row: field hv[N] of algorithm    */
  426.  
  427.     if (!RbfAllocMatrix(2, n, &help))
  428.     {
  429.         ErrMess("RbfInvMatrix: impossible to allocate helpmatrix.\n");
  430.         return KRERR_INSUFFICIENT_MEM;
  431.     }
  432.  
  433.     /* initialize permutation field:                */
  434.  
  435.     for (j = 0; j < n; j++)
  436.         RbfMatrixSetValue(&help, 0, j, (float) j);
  437.  
  438.     /* invert matrix:                        */
  439.     /* notation in comments: m[r,c] means element at row r column c */
  440.  
  441.     for (j = 0; j < n; j++)
  442.     {
  443.         /* looking for pivot:                    */
  444.         /* find row r, where r >= j and m[r,j] is maximal        */
  445.  
  446.         max = fabs(RbfMatrixGetValue(m, j, j));
  447.         r = j;
  448.         for (i = j+1; i < n; i++)
  449.         {
  450.         hi = fabs(RbfMatrixGetValue(m, j, i));
  451.         if (hi > max)
  452.         {
  453.             max = hi;
  454.             r = i;
  455.         }
  456.         }
  457.  
  458.         /* test if matrix is singulary:                */
  459.         /* if true, then return directly and report errorcode    */
  460.  
  461.         if (max == 0.0)
  462.         {
  463.         RbfFreeMatrix(&help);
  464.         return 0;
  465.         }
  466.  
  467.         /* if r != j (r > j) then switch row r with row j and mark    */
  468.         /* this in helpmatrix:                    */
  469.  
  470.         if (r > j)
  471.         {
  472.         for (k = 0; k < n; k++)
  473.         {
  474.             hr = RbfMatrixGetValue(m, j, k);
  475.             RbfMatrixSetValue(m, j, k, RbfMatrixGetValue(m, r, k));
  476.             RbfMatrixSetValue(m, r, k, hr);
  477.         }
  478.         hi = RbfMatrixGetValue(&help, 0, j);
  479.         RbfMatrixSetValue(&help, 0, j, RbfMatrixGetValue(&help, 0, r));
  480.         RbfMatrixSetValue(&help, 0, r, hi);
  481.         }
  482.  
  483.         /* do the transformation:                    */
  484.  
  485.         hr = 1.0 / RbfMatrixGetValue(m, j, j);
  486.         for (i = 0; i < n; i++)
  487.         {
  488.         RbfMatrixSetValue(m, i, j, hr * RbfMatrixGetValue(m, i, j));
  489.         }
  490.         RbfMatrixSetValue(m, j, j, hr);
  491.         for (k = 0; k < n; k++)
  492.         if (k != j)
  493.         {
  494.             for (i = 0; i < n; i++)
  495.             if (i != j)
  496.             {
  497.                 RbfMatrixSetValue(m, i, k,
  498.                     RbfMatrixGetValue(m, i, k) -
  499.                     RbfMatrixGetValue(m, i, j) *
  500.                     RbfMatrixGetValue(m, j, k));
  501.             }
  502.             RbfMatrixSetValue(m, j, k, 
  503.             -hr * RbfMatrixGetValue(m, j, k));
  504.         }
  505.     }
  506.  
  507.     /* now switch back the columns:                    */
  508.  
  509.     for (i = 0; i < n; i++)
  510.     {
  511.         for (k = 0; k < n; k++)
  512.         RbfMatrixSetValue(&help, 1, 
  513.             (int) RbfMatrixGetValue(&help, 0, k),
  514.             RbfMatrixGetValue(m, i, k));
  515.         for (k = 0; k < n; k++)
  516.         RbfMatrixSetValue(m, i, k,
  517.             RbfMatrixGetValue(&help, 1, k));
  518.     }
  519.  
  520.     /* report success:                        */
  521.  
  522.     RbfFreeMatrix(&help);
  523.     return 1;
  524. }
  525.  
  526. #endif
  527.  
  528. /************************************************************************/
  529. /* function RbfMulTranspMatrix:                        */
  530. /* set m1 to m2*m2T. number of rows of m2 must be equal to number of    */
  531. /* rows and columns of m1.                                              */
  532. /************************************************************************/
  533.  
  534. void    RbfMulTranspMatrix(RbfFloatMatrix *m1, RbfFloatMatrix *m2)
  535. {
  536.     int    dest_c;
  537.     int    dest_r;
  538.     int    count;
  539.     register float scalar_product;
  540.  
  541. #ifdef    DEBUG_MODE
  542.     if (m2 -> rows != m1 -> rows ||
  543.         m2 -> rows != m1 -> columns)
  544.     {
  545.         ErrMess("RbfMulTranspMatrix: incompatible matrixes.\n");
  546.     }
  547. #endif
  548.  
  549.     /* notice that m2*m2T is a symmetric matrix!                        */
  550.  
  551.     for (dest_r = 0; dest_r < m1 -> rows; dest_r++)
  552.     {
  553.     for (dest_c = dest_r; dest_c < m1 -> rows; dest_c++)
  554.     {
  555.         scalar_product = 0.0;
  556.         for (count = 0; count < m2 -> columns; count++)
  557.         {
  558.         scalar_product += RbfMatrixGetValue(m2, dest_r, count) *
  559.                       RbfMatrixGetValue(m2, dest_c, count);
  560.         }
  561.         RbfMatrixSetValue(m1, dest_r, dest_c, scalar_product);
  562.         if (dest_r != dest_c)
  563.         RbfMatrixSetValue(m1, dest_c, dest_r, scalar_product);
  564.     }
  565.     }
  566. }
  567.  
  568. /************************************************************************/
  569. /* function RbfMulMatrix:                        */
  570. /* set m1 to m2*m3. number of columns of m2 must be equal to number of  */
  571. /* rows of m3. m1 must be allocated with r = number of rows of m2 and   */
  572. /* c = number of columns of m3.                        */
  573. /************************************************************************/
  574.  
  575. void    RbfMulMatrix(RbfFloatMatrix *m1, RbfFloatMatrix *m2, RbfFloatMatrix *m3)
  576. {
  577.     int    dest_c;
  578.     int    dest_r;
  579.     int    count;
  580.  
  581. #ifdef    DEBUG_MODE
  582.     if (m1 -> rows != m2 -> rows ||
  583.         m1 -> columns != m3 -> columns ||
  584.         m2 -> columns != m3 -> rows)
  585.     {
  586.         ErrMess("RbfMulMatrix: incompatible matrixes.\n");
  587.     }
  588. #endif
  589.  
  590.     /* This seems to be a strange way to multiply two matrices but  */
  591.     /* it prevents the swapper from trashing:                       */
  592.  
  593.     RbfClearMatrix(m1, 0.0);
  594.  
  595.     for (dest_r = 0; dest_r < m1 -> rows; dest_r++)
  596.     {
  597.         for (count = 0; count < m2 -> columns; count++)
  598.         {
  599.         for (dest_c = 0; dest_c < m1 -> columns; dest_c++)
  600.             RbfMatrixSetValue(m1, dest_r, dest_c,
  601.                 RbfMatrixGetValue(m2, dest_r, count) *
  602.                 RbfMatrixGetValue(m3, count, dest_c) +
  603.                 RbfMatrixGetValue(m1, dest_r, dest_c));
  604.         }
  605.     }
  606. }
  607.  
  608. /************************************************************************/
  609. /* set m1 to m2+m3. number of columns of m1, m2 and m3 must be equal.    */
  610. /* number of rows of m1, m2 and m3 must be equal.            */
  611. /************************************************************************/
  612.  
  613. void    RbfAddMatrix(RbfFloatMatrix *m1, RbfFloatMatrix *m2, RbfFloatMatrix *m3)
  614. {
  615.     int    dest_c;
  616.     int    dest_r;
  617.  
  618. #ifdef    DEBUG_MODE
  619.     if (!(m1 -> rows == m2 -> rows &&
  620.           m2 -> rows == m3 -> rows) ||
  621.         !(m1 -> columns == m2 -> columns &&
  622.           m2 -> columns == m3 -> columns)
  623.        )
  624.     {
  625.         ErrMess("RbfAddMatrix: incompatible matrixes.\n");
  626.     }
  627. #endif
  628.  
  629.     for (dest_r = 0; dest_r < m1 -> rows; dest_r++)
  630.         for (dest_c = 0; dest_c < m1 -> rows; dest_c++)
  631.         RbfMatrixSetValue(m1, dest_r, dest_c,
  632.             RbfMatrixGetValue(m2, dest_r, dest_c) +
  633.             RbfMatrixGetValue(m3, dest_r, dest_c));
  634. }
  635.  
  636. /************************************************************************/
  637. /* function RbfPrintMatrix:                        */
  638. /* print out matrix m to stream (file) s                */
  639. /************************************************************************/
  640.  
  641. void    RbfPrintMatrix(RbfFloatMatrix *m, FILE *s)
  642. {
  643.     int    r;
  644.     int    c;
  645.  
  646.     for (r = 0; r < m -> rows; r++)
  647.     {
  648.         for (c = 0; c < m -> columns; c++)
  649.             fprintf(s, "%.4e ",RbfMatrixGetValue(m, r, c));
  650.         fprintf(s,"\n");
  651.     }
  652. }
  653.  
  654. /************************************************************************/
  655. /* function ErrMess:                                                    */
  656. /* print message to stderr                                    */
  657. /************************************************************************/
  658.  
  659. void    ErrMess(char *message)
  660. {
  661.         fprintf(stderr, "%s", message);
  662. }
  663.  
  664. /************************************************************************/
  665. /* end of file                                */
  666. /************************************************************************/
  667.  
  668.  
  669.  
  670.  
  671.  
  672.  
  673.  
  674.  
  675.  
  676.  
  677.  
  678.  
  679.  
  680.  
  681.  
  682.  
  683.